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1. Introduction 

Mobile phone communication is now a key technology across large swathes of the 
world. One of the essential technological ingredients for its success is the ability for 
multiple users to share a single channel (i.e. many people can use the same channel 
to communicate between their mobile phones and a particular base station). Code- 
division multiple access pQ El IS] is a protocol that allows this multiple access to a 
single channel through each user modulating their signal (via so called spreading codes) 
before transmitting to the base station. The base station receives a mixture of these 
modulated signals, combined with aditional channel noise, and the task is then to 
use knowledge of the spreading codes and received signal to reconstruct the original 
information. CDMA has been the subject of several studies in the last few years that 
have utilised the relationship between a model of the communication process and the 
statistical mechanics of fully connected disordered Ising spin systems to examine the 
posterior distribution of the orignal signal using Bayesian inference jU El U\- This 
gives access to maximum a posteriori (MAP) and maximum posterior marginal (MPM) 
decoding. Progress has also been made in terms of algorithms for decoding using message 
passing procedures [HI El EDI- Recently, both density evolution and generating functional 
analysis have been used to analysze the dynamics of some detection algorithms for the 
parallel inference canceller [XT] IT2] . The first of these techniques makes the relatively 
strong approximation of a Gaussian local field, and ignores the Onsager reaction term 
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in the local field, however it is known to generally give relatively good results when 
the detection dynamics converge [TI]. In contrast, the generating functional approach 
[T3| H2] is exact, but the complexity both analytically and numerically increases rapidly 
with the number of time steps considered and thus it is only practically useful for 
examining the first few steps of the dynamics. In the current paper we exploit the 
alternative approach of dynamical replica theory [Til ITSl ITo] . This allows us to treat 
the dynamics of a sequential update detection algorithm (namely the Gibbs sampler) 
working in continuous time, with an analytic approximation scheme that we expect 
to be superior to density evolution and with a numerical effort that increases only 
linearly in time. It was noted in [1] that there exists a spinodal for both MPM and 
MAP decoding, past which the decoding problem has two locally stable solutions, one 
with good performance and one with relatively poor performance. This coexistence has 
practical implications since local search algorithms starting from an initial state with a 
relatively high error rate are closer to the poor solution (at least in the sense of the error 
rate). To go beyond this qualitative argument, however, one really requires dynamical 
tools since the basins of attraction are dynamically defined concepts. The theory which 
we develop in the current paper allows us to examine these concepts. We examine 
the theory for the Gibbs sampler, as a prototype local search algorithm, not because we 
believe that it is necessarily the optimal algorithm for detection and decoding in CDMA 
type problems. The dynamical replica approach we describe here is an approximation, 
but its justification is given by comparison with Monte Carlo simulations. Finally we 
note that the theory we develop here is also applicable with minor modifications to the 
linear Ising perceptron JU] 

2. Model definitions and order parameter evolution equations 

We consider the demodulation problem for the iV-user direct-sequence binary phase- 
shift-keying (DS/BPSK) CDMA system, with the simplifying assumptions that the 
channel noise is additive white Gaussian, chip and symbol timing are perfectly 
synchronized across users and the output power of the users is perfectly equalized by 
power control. For details of the equilibrium statistical mechanical analysis of this model 
please see jU |H| ■ 

We consider N users sending information bits of G {—1,1} Vz = 1, . . . , N. Each 
user i has a binary spreading code {rjj : t = l,...p}, rjj G { — 1,1} so that in 
symbol interval t, user i transmits r^of . We model the speading code sequences to 
be independent quenched random variables with Probfr/' = ±1] = | and take the zero 
mean additive white Gaussian noise \v l : t — 1, . . . ,p} to have variance N/ f3 s . Thus, at 
each chip time step t G {1, . . . ,p} the received signal at the base station y 1 is given by 
jv 



Bayesian inference shows that posterior distribution of the original signal, given the noisy 




(1) 



i=i 
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signal is given by a Gibbs-Boltzmann distribution with temperature (3 and Hamiltonian 
1 N 1 1 p 

h{<t) = -J2 Jaws - E M J * = nJ2 fai * = # E w ( 2 ) 

ij i=l t t=l 

where the signal {y f } and the spreading codes {77*} constitute quenched disorder. The 
temperature (3 is a free control parameter, MAP decoding corresponds to the limit 
(3 — > 00 while MPM decoding corresponds to (3 = (3 S (the Nishimori temperature [B]), 
although in general /3 S may not be known. 

We examine the detection problem by using a spin system with Glauber dynamics to 
model the posterior distribution. To study the dynamical evolution of this distribution 
analytically, we use the techniques of dynamical replica theory [TU HH1 HE] , at the level 
of a three parameter approximation. With Glauber dynamics the time evolution of the 
microscopic state probability distribution Pt(<x) is given by the master equation 

d N 

— p t {cr) = ^\p t (F k cr)w k (F k (T) - p t (cr)w k (cr)] (3) 

fe=i 

with the spin-flip operator F k Q(cr) = $(0"!, . . . , — o~ k , . . . , cr^r) and transition rates w k {cr) 
given in terms of the local alignment fields h k (cr) as 

w k{v) = \\ l ~ o- k tanh(f3h k (er))] h k (cr) = f k - ^ J kj aj (4) 

Conventional demodulation [S] corresponds to taking a® = sgn(f k ) where &° is our 
estimator for cr, the true signal. To improve upon this we take into account correlations 
induced by the spreading code. The dynamics (JHJ) lead asymptotically to the required 
posterior distribution in the high t limit (i.e. the Bayesian posterior distribution). The 
primary performance measure for any demodulator is given by the overlap M between 
the signal er° and the estimate of the signal cr, defined by 

i 

We also use as macroscopic order parameters the internal energy 

E ^ = ? °" ,ij(T ' = f + ^ (Ti ' ,ij ° J ^ 

ij i^j 

(note the similarity to the order parameter r(cr) from ^1) and the contribution due to 
the external fields 

i 

The first order parameter is our performance measure while the latter two give the 
energy from the Hamiltonian. We have chosen to split the energetic term into two 
pieces since we will find that under our assumptions both E and F will evolve according 
to odes containing a relaxation term and a complicated force term. Since, E is quadratic 
in the spins, and F is linear in the spins, if we took E — F (i.e. the energy of the system) 
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to be a single order paramter, then its evolution would follow from the difference of two 
complicated force terms. We found that the degree of analytic complexity was the same 
for either choice. However, by splitting the energy into three terms our approximation 
to Pt(cr) has an extra degree of freedom and thus the approximation is better. Hence, 
at no extra analytic or numerical cost (compared to the standard two order parameter 
theory) we obtain a better approximation to the dynamics. 

Following w e may derive a Kramers- Moyal expansion for the probability 

density V t {M, E, F) = ^ Pt {<r)5[M - M{a)}6[E - E{a)\8[F - F(<r)}, 



dt 



V t {M,E,F) 



V t {M,E,F) 




^a°tanh[/^(<r)] 

i 

l^/4° c (<x)tanh[/^(<r)] 

i 

1 53 /,(«•) tanh^(o-)] 



M 



M.E.Fvt 



+ a-2E 



M.E.Fvt 



— F 



M.E.Fvt 



where we define h 



loci 



(8) 

= Ylj^i Jij a j- I* 1 the thermodynamic limit, on finite timescales, 
only the Liouville term survives in this equation, so that the order parameter triple 
(M, E, F) evolves deterministically according to 



dt 



-M 



-M+/^X)tT°tanh[^( CT )]\ 



M,E,F:t 



dt 



E = -2E + a+ (jjYl /l i° C (°') tanh[/3/i;(<r) 



M.E.Fvt 



M,E,F;t 



where 

</(*)> 



M.E.Fvt 



Her PM5[M - M{<t)}5[E - E{a)}5[F - F{a)}f{c 
Y,(rPt{<T)8[M - M(ct)]S[E - E(tr)]5[F - F{*)] 



(9) 
(10) 

(11) 
(12) 



This flow equation is still exact in the thermodynamic limit. However, to move to a 
practical representation, i.e. one that does not depend on the microstate probability 
distribution pt(er), we make the assumptions underlying dynamical replica theory 
|14[ ITT)] : that the observables are self-averaging with respect to the realisation of the 
disorder and initial conditions and that we may appoximate the microscopic probability 
distribution by the maximum entropy distribution given the values of our observables. 
Thus, with this approximation the microstate probability drops out of our equations 
and we may then use the replica technique to remove the unpleasant fraction in ([12)1 . 
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via 



V ^ (T\...,(T n a=l 

which we proceed to do in the following section. 
3. Replica calculation of the flow 

Using site equivalence under the disorder average, the objects we need to calculate can 
be expressed as 



D(h) = (a^6[h-h 1 (a-)]) M ,E,F (14) 
D(f, h) = (S[f - h]5[h - h[° c {<r)}) M ^ F (15) 

From the definition h\(cr) = fx — h l ° c (cr) it is trivial to obtain the first distribution from 
the second, so we focus on the calculation of the second term in the current section. We 
introduce Fourier representations for the delta functions over (M, E, F) constraining 
the distribution in (fTH|) giving conjugate parameters {M a , E a , F a }. To average over 
the disordered signal y t , we generate its correct measure by using the partition function 
with the true value of the noise (3 S j^J. We also write the delta functions over the fields 
in Fourer representation 

dhdf 



sif-^yrfMh-hx^)]^ J 



^ + vf - % E - 1 E E ^» \ (i6) 

t=i j>i t=i 



(2tt) 2 

exp 

To get the correct scaling, we change variables r* = y*/ \J N and then introduce 

1 ~ n n / dvtaddta ex p { ' wta%ta - 7^ to e 4°< \ (17) 

t=l cs=0 L »>1 J 

We then find 

D(f, h) ~ lim / ^ J][dM a dE a dF a ] JJdr* J][dt; to d{) to ]e i ^ +i//+iiV S^ MA:fa+ ^ Q ) 

~^ ^ ' a t=l ta 

e-'E^o^Ei^ ^e-^ELi^C/V+^^-^E^^EiM^f 

<T°,...,<T™ ?7 

e "^7 Eat ^M'f E t uMCrW )-^ E«77*af «*» 

By ~ we mean that this is correct up to irrelevant normalisation constants which can 
be recovered later by dividing through by (1)m,e,f or requring f dfdhD(f,h) = 1. 
Performing the trace over rj in the last line, and then expanding the resultant formulas, 
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exp J — ^[/ViV - v t0 ) - i o-i(E a v ta + #V) - ifr* - ihv n f 

L t a 



- E 



2iV 

i,i>l,a,/3 



We then introduce the Edwards-Anderson order parameter: 
1 = / ]^d^ T d^ r e iE ^ E *>i<<- WE ^ r ^ 



(19) 



leading to 



D(f, h) ~ lim / ^ ]7[dM a d£ a dF a ] fftdvdg^le^^ (20) 

CT, 

JJ J y drdvG(r, v) exp ^[/3 s a°(r - v°) - i ^ a"(E a v a + F a r) - ifr - ihv 1 ] 



t 

Where 



Qfa v ) = e -^'^-^(r-v ) 2 -^ £ a ^ a ^-iE Q ^21) 

It is instructive to view (|20|) as a measure over (f,h), with free parameters which are 
subsequently averaged over a saddle point measure. Thus, the free parameters within 
the order one measure take their saddle point values. The saddle point surface itself is 
given by 

V = i ]T(MM Q + EE a + FF a ) + log eiEQ0 " [ E r «^^ ~ 1 E 

a CT a/3 a>0 

+a log y drdve-^ lv -^-^ 2 -^ ^ t «" i ^° ^ (22) 
4. Replica symmetric saddle points 

We make the replica symmetric assumptions, r 0r = r, r pr = R, q 0a = m and q a p = q and 
(M a ,E a ,F a ) = (M,E,F) V«. Then, similarly to the original equilibrium calculation 
jU El IE] , we have a saddle point surface 

If 1 

= lim = / £>z log 2 cosh( V^Rz + r) - rm i?(l - g) + M(M - m) + EE 

n^o n J 2 

+FF + odog y drdvexp \ -^vq _1 v - y(r - w ) 2 - 7^ XX ~ ( ( 23 ) 
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We wish evaluate the integrals / in the second line of the above so that we can take the 
limit n — > 0. It is convenient to change variables to 



tm 



v 



U\ll 

v a = z a \J\ - q-ty/q 



(24) 
(25) 



where u, t, {z a } are zero mean, uncorrelated, unit variance Gaussian random variables. 

1 12 

Then using the shorthand Dx = (2n)~2e~2 x we have 

2" 



DrDtDu exp 



ft 
2 



tm 



(26) 



x 



Dz exp 



E 

-^{zy/l^q - t^fq) 2 - Fr(zy/T^ - ty/q) 



Although cumbersome, these integrals are straightforward and give 

^ 2 - 1 



m 



\/ba — 



where 



a = 



Ps 



b = 



i + 



9 

Am 2 



nF 2 (l -q) 
l + E(l-q) 

nEq 



l + ft(l 
/? — 



+ 



nFy/q 



i + A(i- 



(27) 

(28) 
(29) 
(30) 



Giving the replica symmetric saddle point surface 

\fr RS = j Dz log 2 cosh{\/Rz + r)-rm- ^R(l - q) + M(M -m) + EE + FF 

p s {Eq + 2mF - F 2 (l - q)) - F 2 (l - q) 



a 
2 



log[l + £(!-<?)] + 



p s [l + E(l-q)} 
Extremizing this surface we find the saddle point equations 



(31) 



r,R: 
M, E 

F : 

m, q : 



m 



J Dzt&nh(\/Rz + r) g = J D z t&nh 2 (VRz + r) 



(32) 



M = m ^ = a J £[1 + £(1 - g) 2 ] - 2mM(l - q) + F 2 (l - q) 2 (l + ft) 



a 



P s m-(l + p s )F(l-q) 



r + M 



(3 s [l + E(l-q)] 

-aF 
~ l + E(l-q) 



R 



(3 s [l + E{l-q)Y 



a[F 2 (l + ft' 1 ) + E(Eq + 2mF)} 
" [l + E{l-q)Y 



(33) 
(34) 
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Equilibrium corresponds to M = and E = —F = j3 since in this case the magnetisation 
is unconstrained while the variable conjugate to the energy is, of course, the temperature. 
It is a useful check to see that in this case the saddle point equations for r, R, m and q 
resort to their standard equilibrium expressions (as then our maximum entropy measure 
is just the equilibrium measure at temperature (3). 

4-1. Reduction of the saddle point equations 

We have to express the values of the unknown conjugate order parameters given the 
values of the true order parameters. By eliminating F from the saddle point equations 
for F and E we obtain a quadratic equation in E for which the physical solution (utilising 
the fact that in equilibrium we require E — 0) is given by 

E{q) = ^ )2 { " M 1 - q) + A(l + A)(l " V?] (35) 
-vW -q)+ Ps(l + Ps)(l ~ q) 2 } 2 ~ 4rf(l - q) 2 [d + (3 S {1 + (3 S ) - M 2 & 2 ]} 
d 



F 2 P 2 2(3 s (l + f3 s )E 



a 2 a 



Insertion of (|33|) into (|3*3|) gives us F(q) which in turn implies R(q). To obtain r(q) we 
have to use the implicit equation for r 



M = J Dzt&nh(^/R(q)z + 



(36) 



We have resolved all free parameters into bare functions of q and hence have the relatively 
straightforward one dimensional problem 



q 



J Dztswh 2 (y/R(q)z + r(q)) (37) 



which we solve numerically. 

5. Derivation of the force terms in the order parameter flow 

We are now in a position to calculate D(h) and D(f, h). We abuse notation slightly by 
using M, E, F, q when we mean their values taken in the saddle point given the values 
of M, E and F rather than seeing them as variables. With that taken into account, we 
may write 

D(h) = hm f ^e lhh Y aoe-M^o^+m*) ( 38 ) 
n-*o / 2n ^— ' 

D(f,h) = hm I f^-^h+ifiy e -ME Q>0 ^+fA(o-) (39) 
n^o J (2n) 2 
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A(< 



/ drdvG(r, v) (3 s (r - v°)a° - J2 a ° a (Fr + Ev a ) - ihv 1 - ifr 



1 2 



J drdvG(r, v) 



(41) 
(42) 



5.1. Evolution of the magnetisation 



We concentrate first on S. It is convenient to introduce the shorthand for averages over 
the measure 



(...) 



G 



J drdvG(r, v) 



(43) 



J drdvG(r, v) 
with which we define the shorthand: 

3io = aP.((r - v°)(r - ^)> G g 0a = a(3 s ((r - v°)(Fr + Ev a )) G + M (44) 

g n = «((r - v l ) 2 ) G g a p = a((Fr + Ev a )(Fr + Ev p )) G (45) 

g la = a((r-v 1 )(Fr + Ev a )) G Va > 1 #m = a((r — v 1 ) (Fr + Ev^q (46) 

We calculate these factors in Appendix A[ for now we merely note that g a p = R while 
—9oa = t. It is possible to ignore constants such as /3 2 ((r — v ) 2 ) G since they may be 
dealt with via overall normalization of the measure D(h). Then 



Din) ~ lim / — e 2 

Since we only require D(h) in a term of the form J dhD(h)tax\h(/3h) it is possible to 
make the gauge transformation x,h,h — > a x, a h, cr h to remove the dependence on ct . 
We then perform the trace over the other spins to obtain 



D(h) 



lim / ^ e -V*** 
♦0 ./ 2vr 



Dxe~ ih9l ° cosh[ihg in + x\[R~ + r] cosh n-:l [i^ la + x\[R- 



In order to move our integration contour so that we can perform these integrals we 
expand the first cosh function, 

cosh[i/i(7iii + xVR + r] (47) 
= cosh(iftA) cosh[i/igi a + xy/~R + r] + sinh(iftA) sinh[i/igi Q + xVR + r] 

where A = gm — g\ a . The integral over the contribution to D(h) containing only cosh 
functions is relatively straightforward and gives 



1 



exp 



(h-gxo + Af 



2g 



11 



+ exp 



(h-g 10 -Af 



'2g 



11 



(48) 



For the term containing the sinh functions from (j47)l . we shift the integration variable 



x — > x — ih^= and take the limit n —>■ to obtain, 



2vr 



sinh(iftA) / £)xe lx?i ^ tanh(xv^R + r) 
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We proceed by writing sinh(iftA) = |(e 1,lA — e lhA ) and we first treat the term with 
containing e lhA , integrating over h: 



47T ' 



2J27ra(g n - 9 -j* 



Dxexp 



2 ^ii ~ t 



-{h- g 10 + A + x 



9la n 2 



R 



(49) 

tanh(xv^R + r) 



By a careful change of variables this reduces to 



1 fe 2 

= e 29 n / .Dxtanh 



rd R 9n-9L _9_ 1 a h + r 
9n 9n 



VZngii 

The case with e~ lAh follows identically, so finally we find that 
D(h) 



(50) 



e -2Hl(^-Sio+A) 



V^9 



11 



Dx tanh 



x 1 



(Rgn - gl 



9u 



9la 

9n 



{h-g w + A)+r 



e-^C'-flio-A) 



V2Tg- 



ii 



tanh 



( Rgu - gjc 

9n 



9la 

9n 



(h - g 10 - A) + r 



(51) 



We wish to consider the term J dhD(h) tanh(/3/i) which is required to calculate the force 
term in the differential equation ®. This can be most easily effected by first making 
the transformation h —>■ h + g±o =F A as required, followed by h — > J~g~x\h. We then have 



dhD(h) tanh(/3/i) 



DhDx < 1 + tanh 



+ - / DhDx < 1 - tanh 



(52) 



R9ll ~ 9la 9lah 



9n 



>9ii 



tanh[^(y^/i + 5-10 - A)] 



, ^11 - 9la 9lah 

x\ I h r 



#11 



'011 



tanh[/3(^iT/i + £io + A)] 



(53) 



We then rotate the Gaussian integration variables which leads to our final result 

&hD{h) tanh(/%) (54) 

1 



DuDv { 1 + tanh 



Ru + r 



tanh[/3( 



- ^j=u + g l0 - A)] + 
v R 



~J DhDx^l -tanh y/Ru + r }tanh[/3( 



^11 - glc 

R 



R9ll — 9la 9la a \ 1 

^ u + 0io + A)J 



i? 



5.£. Fixed points for the magnetisation in equilibrium 



As a useful test of our analysis thus far, we show that in equilibrium, the equilibrium 
value of the magnetisation is a fixed point of our dynamics. In equilibrium we have 
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M = and E = — F = (3 as argued above. So we can read off the equilibrium values 
of our factors g... from | Appendix A" In equilibrium —f3g\ a = R and (3gio = r. Further, 



(3 9 ia = —R while Pg\\ = —gin, whence 

PA = P(jg m - g Xa ) = R- f3 2 g n (55) 

and so 

J dhD(h) taoh(/3/i) = J Dutanh[VRu + r] (56) 
according to the identity 

tanh(-u) = -[1 — tanh(u)] J Dytanh(u + yz — z 2 ) 

+ i[l + tanh(u)] y J Dytanh(w + ?/2 + 2 2 ) (57) 

of [T^j. Since our saddle point equations (j32H34|) for m,q,R,r in equilibrium are 
equivalent to the equilibrium saddle point equations we see that in equilibrium, one fixed 
point of our differential equation for the magnetisation Q is given by the equilibrium 
magnetisation. 

5. 3. Evolution of the energetic terms E and F 

We now turn to D(f, h), for which we need the evaluation of A(er) as defined in (|4*T|) . 
On top of our previous definitions we now introduce 

g rr = a(r 2 ) G g vv = a(vf) G (58) 

g rv = a(rvi) G g ra = a(r(Er + Ev a )) G (59) 

g va = a(vi{Fr + Ev a )) G Va > 1 g vl = a{vx{Fr + Ev x ))g (60) 

g 0r = ap s (r(r - v )) G g 0v = af3 s (vi(r - v )) G (61) 

which are also calculated in appendix |Appendix A| Using these g factors we may write 

D(f h) ~ lim ( d ^ d ^ c -^^-A E£ -/hg rM +ifeh+i//* jhW^vl+Eoi ^g va -a og0v ] 
U ' ' n->0 J (27r) 2 ^ 

J Dxe xy ^^ aCra+rcr °^' aCra 'e~ i tt cro9r0 ~^<* cr <*9r a ] ^g2) 
which we may rewrite as 

D(f,h) ~ lim [ ^ e -^-^-fh 9r «+ihh + iff f Dx y e -ia (h 90v+ f 9 or) (63) 

n->0j (2n) 2 J 

cosh[i/t(7„i + ifg ra + xVR + a°r] cosh™^ 1 [ihg va + ifg ra + xVR + cr°r] 
Following a similar procedure to that used in § 15.11 we write 

cosh[i/i^i + ifg ra + x\fR + a°r] = cosh[— ihA] cosh[ihg va + ifg ra + x\f~R + cr°r] (64) 

+ sinh[— ihA] smh[ihg va + ifg ra + xVR + a°r] 
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since g v \ — g va = —A. It is convenient to use the change of variable h — *> h + go v o~o ± aA 
and / — > f + gorO'o. The terms in (jo^|) containing the cosh terms from ()64|) can be 
treated in a similar manner to § 15.11 and contribute 

y 



-g 2(grrgvv—9rv) 2(g rr gvv — g£ v ) (grrgvv— grv) 



_. (65) 

Now we tackle the more tricky sinh terms from (jo^|) which requires the shift of the x 
integral in the complex plane but then gives 



E 



TO 



2(2tt) 5 



e 2 2 

2(2vr) z 

/ — <.</<■: -IT ' — ~7^ ; '" ;f;I /v 



(S- - ^ ) - V (9— *W-)-fh(grv- 



H / Dxtanhfxv^R + i(hg va + fg ra ) + ^o"o] 



2(2tt) 1 



77? 



(66) 



x j Dx tanhfxv^R + ra ] 



Proceeding with the integrals over h, f and after some rearrangement we find, 

] grrft 2 gyyf 2 L grvfh 



E 



cro 



2 

rv 



-g 2(grrgvv—g r v) 2(grrgvv—g r v) (grrgvv—g r v) 



2ny/g rr g vv - g 

X ( Dx tanh[a; I ^ \9rr9va 9ra.grv) ^ (gvvgra gva9rv) j. ^ j 



(9rr9vv g 



C (g>rrgw — g1 v 
Putting this all together we have 

dfdh 

2 



2 ) 

rv i 



(67) 



dfdhD(f, h)fttmh[P(f -h)} = J2 J 



47r/ 



g rr g vv g r 



9rrh 



gw f 



+ - 



grvfh 



g ^(grrgvv—grv) ^(grr gw —grv) (grrgvv— 9rv) 



i / , i I R {grrgva 9ra9rv) , \9w9ra 9va9i -v j ,• i 

1 + <Ti / Dx tanh[x\/ — t^* 1 7 ~ — / + rcr oJ 



(68) 



C {grrgvv 9rv) {grrgvv 9rv) 

(/ + o- g 0r ) tanh[^(/ + g 0r cr - {h + g 0v a + o"iA))] 

Rescaling our Gaussian variables in order to write the measures as standard Gaussian 
measures we finally have the force term required for ([TTf 



dfdhD{f,h)ftanh[P{f-h)\ (69) 
= ^J2j DyDzfi + nJ Dxt&nh{A + B)]Ct&nh{f3{C - D)] 



with 



.4 



x i 



R{9rr9vv 9rv) 9va9rr 9ra9vv ^9ra9va9r 



B = ro~Q — a/^ 1 

C = JgVrV + gorO-o 



{9rrgw g"^v) 

grrgva QraQrv 



\fg rr gw - g 



+ g ra y 



(70) 

(71) 
(72) 
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updates/spin 

Figure 1. We compare the results for the magnetisation M of solving our order 
parameter equations numerically (solid line) with Monte Carlo simulations (dotted 
line) for MPM decoding at f3 s — (3 — 2 and a = 2. The Monte Carlo simulations are 
performed with system size 2 000 averaged over 50 samples (the standard deviation is 
less than 10~ 3 ). 



D 



(9rr9 



9rv) I 9rv \ 

9r 



while with these definitions, the force term required for (|T0|) is 
dfdhD(f, h)htanh[(3(f - h)) 

= ~ J DyDz[l + a 1 J Dx tanh(A + B)\D tanh[/3(C - D)\ 



(73) 



(74) 



6. Order parameter flow and comparison with simulations 

We now have a closed set of equations describing deterministic order parameter flow 
They are 

d 



dt 



M = -M + / dhD(h) tanh((3h) 



= -2E + a + J dfdhD(f, h)ht&nh[(3(f - h)} 
±F = -F + J dfdhD(f,h)fUnh[P(f-h)} 



(75) 
(76) 
(77) 



where D(h) and D(f,h) at any instant depend on the triple (M,E,F) (as well as the 
statistics of the quenched disorder via (3 S and a). 

As an initial basic test of our theory, we compare the flow described by the solution 
of our order parameter equations to Monte Carlo simulations of the original Glauber 
dynamics in figure H The temperature is relatively high compared to realistic values 
but we see that at least in this regime we have excellent agreement between the theory 
and simulations, justifying our assumptions and validating our method. 



Dynamical replica theoretic analysis of CDMA detection dynamics 



14 



6.1. Phase coexistence and basins of attractions 

It was noted by Tanaka 01 IHj that for certain parameter regimes a spinodal would be 
encountered. This has drastic implications for the decoding problem, since if there are 
metastable states with a high bit error rate, local search algorithms will fail to find 
the low bit error rate solution. To go beyond these arguments a dynamic approach is 
required and we have now provided one. It is possible to look at flow in parameter space 
and see explicitly the basins of attraction for both the good and bad solutions. Spinodals 
occur for both MAP decoding and MPM decoding. We focus on the latter since the 
MPM decoding is optimal in terms of the bit error rate. Since it is at Nishimori's 
temperature, for equilibrium states at least, there are no complications due to replica 
symmetry breaking (we cannot guarantee this for the dynamical saddle point at present, 
but we hope to investigate this further in a later work). 

There are a variety of initial conditions with which we could start, however two are 
of particular note. The first is the random initial state (i.e. cxj(0) = ±1 with probability 
|). It is straightforward to derive the initial states of our order parameters as 

M RIS = F RIS = E RIS = ^ (78) 

A second important initial state is that given by the conventional demodulator |3| , with 
<7j(0) = sgn(/j). We derive the values of our order parameters in appendix |Appendix B[ 
which are 



M CD = Erf 



a 



2(1 + 



CD 



Dz 



a + y/a(l + Pi 1 )* (79) 



E C d = | [1 + 2M CDX + X 2 (l + P; 2 )} X = \j na ^ + p-if ( 8 °) 

where Erf(x) = ^= J^dte - * 2 . 

In figure [21 we plot the flow of the dynamics in the energy - overlap plane. We see 
that as expected the random initial state (RIS) and conventional demodulator (CD) 
state both flow into the poor attractor. The agreement between theory and simulations 
is quite reasonable up to some low energy point where we believe replica symmetry 
breaking (RSB) starts to play an important role in the simulations (although our current 
replica symmetric theory is unable to take it into account). We also note that starting 
from the CD state the overlap gets worse before improving to a state better than the CD 
in the theoretical approach, while in the simulations, it never recovers its initial overlap. 
Thus, starting from CD and running a local search algorithm can significantly decrease 
performance levels in a practical setting. We also see that the basin of attraction for the 
good final state is significantly smaller than that for the poor final state, the boundary 
is well over half way between the two states in (E — F, M) space. This also shows how 
distance in (E — F, M) space is an unreliable guide to the direction of flow towards 
attractors. 
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Figure 2. We examine the flow through phase space projected onto the energy (E-F) 
- overlap (M) plane. We both solve our order parameter equations numerically (solid 
line) and compare against Monte Carlo simulations (dotted line) for the two important 
initial staes CD and RIS. We work with MPM decoding at /3 S = = 17 and a = 0.5. 
The crosses mark the two stable solutions to the equilibrium problem. The Monte 
Carlo simulations are performed with system size 2 000 averaged over 50 samples. We 
ran both simulations and our theory for 1 000 updates per spin. We see that for 
this parameter regime there is phase coexistence. The line labelled RIS starts from a 
random initial state, while the one labelled CD has the conventional demodulator as 
its initial state. Both flow into the poor attractor. We also show flow starting from 
some other (non-physical) initial states to get a better idea of the basins of attraction 
for the two phases. Note that all flow is from top to bottom as the dynamics lowers 
the energy of the system. 



7. Conclusions 

CDMA is an important standard used in modern mobile communications. Tools 
from statistical physics have provided and will continue to provide useful ways of 
examining the detection problem. Here we have developed and used dynamical replica 
theory to study the dynamics of the detection problem for a prototypical local search 
algorithm, namely the Gibbs sampler under Glauber dynamics. Although this approach 
is only an analytic approximation, it provides a useful counterpart to both density 
evolution and generating functional analysis as a tool for examining dynamic rather than 
equilibrium properties. As we have seen in comparison with Monte Carlo simulations the 
approximation is a reasonably good one. We have also calculated the basin of attraction 
for a particular set of parameters, in the region where there is phase coexistence. As 
expected we have seen that the practically available initial states, CD and RIS, both 
flow to the poor solution. We have also seen that in this case the overlap decreases 
initially from the CD state with our search algorithm and that the basin of attraction 
for the good solution is relatively small in energy-overlap space. One obvious extension 
of this work would be to increase the order parameter set to improve the level of 
approximation. Since the number of updates per spin required to visit interesting 
regimes is of the order of 10 3 , this would have to be done in a way that was compatible 
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with practical numerical solution. Another interesting problem is examining the role 
of replica symmetry breaking on the dynamics. Finally, a project that we are already 
working on is using the dynamical replica approach for parallel update dynamics where 
we can compare its predictions to the exact theory of generating functional analysis |12j . 

Appendix A. Calculation of factors defining D[h) and D(f,h) 

The various factors defined in (j44H4(i| are all simple combinations of 

(r% (rv°) G (rv a ) G (v°v a ) G (v a v?) G (v a v a ) G (A.l) 

which are moments of the measure G(r, v) (|4"2jl The calculation is not dissimilar to that 
carried out in the equilibrium calcultion; we have to just be careful with the algebra 
over various Gaussian integrals. We find, 

(v g v a )t 1 ~ g , <? + (1 ~ g) 2 ^ 1 + ft" 1 ) - 2(1 ~ a)Fm 

l + [l+i(l-g)] 2 

(r 2 h =(l + /3; 1 ) (A.3) 

(rv°) G = 1 (A.4) 

(r V ") G = '"" il • ^ (A.5) 

[l + E(l-q)\ 

(v°v») G = m -(}-«) P (A.6) 
[l+E{l-q)\ 

WV = ^ + (I -qy^ + PT 1 ) -2(1 -q)Fm 

[l+E(l- q )]z 

Hence the various factors are simply 

l + (l-q)(E + F) _ c . 

g w = a g 0a = r g a0 = H (A. 8) 

[l+E(l-q)} P 

and 

(3- 1 + 2 - 2m + E(l - qf + 2(1 - q)(E + F)((3~ l + 1 - m) + (fc 1 + 1)(1 - g) 2 (£ + F) 2 
On = a * 

[l + S(l-g)] 2 

We then have 

(1 + fr)F + (E- F)m -Eq + (E + F)(l- q)(Em + F(l + fc 1 )) 

9ia = ct (A. 9) 

[l + E(l-q)} 2 

which is — R in equilibrium. Then 

A ^-*E0-q L F m + Eq-(l-q) [m EF + F 2 (l + ^) ] 

l + E(l-q) [l + E(l-q)} 2 

and 

F(l + + Em n .. . . . /*s 

o rQ = a — 9ro = a(J s {{r ) G - {rv ) G ) = a (A.ll) 

1 + £(1 -g) 

fto = otf3 s ({rv x ) G - (viv a ) G ) (A. 12) 
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Some of these factors only appear in certain combinations, however, although we 
have tried we have not made any significant simplification through further algebraic 
manipulation. 

Appendix B. Calculation of our order parameters for the conventional 
demodulator 



In this appendix, to lighten notation slightly we assume that the original message was 
cr while our estimator is <x. The conventional demodulator sets G{ = sgn(/j) with 

= ^ £ Vtvfr + ± £ r,y (B.l) 

Since v l ~ J\f(0,N/(3 s ) and we assume that {77*, a^v*} are all mutually independent 
random variables, we can use the central limit theorem to treat the second and third 
terms in (jB.ljl . So the error probability for a single bit is the probability that a Gaussian 
7V(0, 1) random variable is less than — a/«/(1 + which gives 



Mi 



CD 



Erf 



a 



2(l + /3 7 i)J 



We can also see that F CD = fi s S n (fi) =^\fi\ is g iven by 



F CD = j Dz 
Now, we can write Eqd as 



E, 



CD 



CD\2 



ni 



CD 



(B.2) 



(B-3) 



(B.4) 



Now, 



bi = ritaga[oaTi + ^ £ ViVfa + £^* 



(B 
(B 



+^ £ (ll fa + A 26 + ^ £ fate + jf £ tiA + °( N 

i \ j / \ t^sjjH t^s / 

where we have thrown away irrelevant terms and expanded the sgn function in a Taylor 
expansion (it is possible, although longer to check we obtain the same result without 
using this trick). We define 



7 



X 



^£25(^+1 Yl rivh + ^Y,viA 



(B.7) 
(B.8) 
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2 



e w+f>7 1 ) (B.9) 



a 



(B.10) 



j 

So up to irrelevant constants, 



mf D = 7 * + k*x + y/Np^xzi 



(B.H) 



where above and in the following z±, z%,Z3 are Af(0, 1) random variables. The slight 
complication is that 7* and are correlated since the latter is a trace over njai while 
the former is a trace over essentially rj\ai - thus Mod determines the degree of correlation 
and we find 
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